Analysis of the elasto-plastic response of a polygonal packing 
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We investigate the constitutive response of two-dimensional packed samples of polygons using 
molecular dynamics simulation. The incremental elasto-plastic response is examined in the pre- 
failure regime. Besides the Young modulus and the Poisson ratio, an additional parameter must be 
included, which takes into account the anisotropy of the elastic response. The plastic deformations 
are described by the introduction of the yield and the flow directions. These directions do not agree, 
which reproduces the non-associated feature of realistic soils. In order to detect the yield surface, 
different loading-unloading-reloading tests were performed. During the reload path, it is found that 
the yielding develops continuously with the amplitude of loading, which does not allow to identify 
a purely elastic regime. 



I. INTRODUCTION 

Traditionally, the quasi-static deformation of soils has 
been described by using constitutive laws. They are em- 
pirical relations between the stress and the strain involv- 
ing a certain number of material parameters, which, in 
the simplest models, can be measured in experimental 
tests However, the more sophisticated models in- 

volve so many parameters that their direct experimental 
meaning and their identification becomes impossible. 

In the last years the numerical simulations have been 
used as an alternative to study the behavior of soils. Usu- 
ally, disks or spheres are used in order to capture the 
granularity of the materials [?J . The simplicity of their 
geometry allows to reduce the computer time of calcula- 
tions. However, they do not take into account the diver- 
sity of shapes of the grains in realistic materials. 

A more detailed description is presented here by using 
randomly generated convex polygons. The interaction 
between the polygons can be handled by letting the poly- 
gons interpenetrate each other and calculating the force 
as a function of their overlap ||. This approach has 
been successfully applied to model different processes, 
like fragmentation , damage Jl2|,[l0| , strain localiza- 
tion and earthquakes The contribution of this work 
is to determine the constitutive relation of this discrete 
model material in the regime of quasi-static deforma- 
tions. The results show that simple mechanical laws at 
the grain level are able to reproduce the complex macro- 
scopic behavior of the deformation of soils. 

The details of the particle model are presented in Sec. 
[n|. In addition to the normal contact force mentioned 
above, the tangential contact force law is implemented 
by a Coulomb friction criterion, and the boundary condi- 
tions are modeled by the introduction of a flexible mem- 



brane that allows to fix the stress value. The calculation 



of the constitutive relations is presented in Sec. III. Wc 



discuss the results in the framework of the classical theory 
of clasto-plasticity. The summary and the perspectives 



of this work are presented in Sec. [V 



II. MODEL 

The polygons representing the particles of this model 
are generated using a simple version of the Voronoi tessel- 
lation: First, we choose a random point in each cell of a 
regular square lattice, then each polygon is constructed 
assigning to each point that part of the plane that is 
nearer to it than to any other point. Each polygon is 
subjected to interparticle contact forces and boundary 
forces as we explain below. 

When two polygons overlap, two contact points appear 
form the intersection of their edges. The contact line is 
defined by the segment connecting these two intersection 
points. The contact force is calculated as 



/ c = k n Ax c n n c 



k t Ax\t c 



(1) 



here n c and t c denotes the normal and tangential unitary 
vectors with respect to the contact line, and k n and k t 
are the stiffness in the respective directions. The over- 
lapping length Ax„ is the ratio between the overlap area 
of the polygons and the length of their contact line. Ax^ 
defines the elastic tangential displacement of the contact, 
that is given by the time integral starting at the begin of 
the contact 



Ax c t 



(2) 



Where is the Heaviside function and v$ denotes the 
tangential component of the relative velocity at the 
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contact. v° depends on the linear velocity Vi and angular 
velocity uji of the particles in contact according to 



Ui x £i + ujj x £j 



(3) 



The branch vector li connects the center of mass of par- 
ticle i with the point of application of the contact force. 
This point is taken as the center of mass of the overlap- 
ping polygon. Eq. (^) defines a limit of elasticity in the 
contact force. When the contact force satisfies / t c = [if£ 
the contact slips, giving rise to a plastic deformation. 

The external forces are applied on the boundary 
through a flexible membrane which surrounds the sam- 
ple. Such a membrane is calculated using an iterative al- 
gorithm, that selects the segments of the external contour 
whose bending angle is smaller than a threshold angle 6 t h 
|p^| . On each selected segment T = Axixi + Ax^x^, we 
apply an external force of the form 



/ = — (JiAx^sti + (J3AX1X3. 



(4) 



Here x\ and £3 are the unit vectors of the Cartesian 
coordinate system. o~\ and 03 are the components of the 
stress we w ant to apply on the sample, as it is presented 
in Sec. fill 




FIG. 1. Schematic plot of the membrane obtained with 
threshold bending angle 7r/2. 

The contact forces and the boundary forces are in- 
serted in Newton's equations of motion which is solved 
numerically using a predictor-corrector algorithm. In or- 
der to enhance the stability of the numerical method and 
to allow for rapid relaxation, some viscous forces are in- 
cluded both in the contacts and in the boundaries: 



fv = -m{~f n v c n h c + ~f t Vti c ), 



f v = -mcitVi 



(5) 



The contribution of these forces is almost negligible in 
the quasi-static regime where velocities are small. They 
are included only to reduce the acoustic waves emitted 



when the system goes from one equilibrium state to the 
other, m = (l/rrii + 1/rrij)^ 1 is the effective mass of the 
particles in contact, and m, is the mass of the particle i 
in contact with the membrane. 

There are three characteristic times in the simulation: 
The relaxation time t r — l/"f n , the loading time <o an d 
the characteristic period of oscillation t s = ^Jk n /mty. 
Here mo is the mean mass of the polygons. This leads 
to a minimum set of dimensionless parameters, whose 
selected values are shown in the table. 



variable 


ratio 


default value 


time of relaxation 


trjt s 


0.1 


time of loading 


to/ts 


1250 


friction coefficient 




0.25 


stiffness ratio 


h/K = Jt/ln 


0.33 


bending angle 


6th 


tt/4 



III. CONTINUOUS RELATIONS 

The characterization of the macroscopic state of a 
granular material in static equilibrium is usually given 
by the Cauchy stress tensor. The derivation of this ten- 
sor over a representative volume [|l4fl leads to 



1 



A^ 

b 



x i f j 



(6) 



The sub-scripts i and j in Eq. (g) denote the com- 
ponents of vectors and tensors. Here x b is the point of 
application of the boundary force f 3 . This force is de- 
fined in Eq. . A is the area enclosed by the boundary. 
The sum goes over all the boundary forces of the sample. 
Inserting Eq. in Eq. (||) leads to 



c = ~r 



-a, J2 b x b Ay b 

-ai E b y"^y b 



c b Ax b 



°3 J2b 



(7) 



These sums can be converted into integrals over closed 
loops. Then, the calculation of such integrals leads to 



o 1 

CT 3 



(8) 



Thus, the stress eigensystem coincides with the Cartesian 
coordinate system used. We can reduce the notation in- 
troducing the pressure p and the shear stress q in the 
components of the stress vector 



(9) 



In the same way, the incremental strain tensor can be 
calculated as the average of the displacement gradient 
over the area of the sample. It has been shown jl5| that 
this average can be transformed into a sum over bound- 
ary segments of the sample 
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dei 



(10) 



Here N is the 90° counterclockwise rotation of the 
boundary segment T. The displacement of the segment 
Ax b is calculated from the linear displacement Ax and 
the angular rotation Ac/) of the polygon, according to 



Ax 6 = Ax + A<p x t 



(11) 



The vector I connects the center of the segment with 
the center of mass of the polygon to which it belongs. 

The eigenvalues dei, de% of the symmetric part of de,j 
define the volumetric and shear component of the strain 
as the components of the incremental strain vector. 



de = 



de v 




dei + de 3 


de-y 




dei — de3 



(12) 



From the macro- mechanic point of view, each state 
of the sample is related to a single point in the stress 
space, and the quasi-static evolution of the system is 
represented by the movement of this point in the stress 
space. The resulting deformation during the transition 
from stress state a to a + da is given by the incremen- 
tal strain de. In advance, let us separate the incremental 
stress in its elastic (recoverable) and plastic (irrecover- 
able) components: 



de = de e + de p . 



(13) 



Following the procedure proposed by Bardet |?| both 
components can be obtained as it is shown in Fig. 3. Ini- 
tially, the sample is in the stress state a. Loading the 
sample from a to a + da the strain increment de is ob- 
tained. Then the sample is unloaded, back to the original 
a, and one finds a remaining strain de p , that corresponds 
to the plastic component of the incremental strain. For 
small stress increments the unloaded stress-strain path is 
almost elastic. Thus, the difference de e — de — de p can 
be taken as the elastic component of the strain. 

This procedure is implemented on one sample choosing 
different stress directions. Fig. || shows the load-unload 
stress paths and the corresponding strain response when 
an initial stress state with q = p/4 is chosen, correspond- 
ing to ai = bp/A and 03 = 3j?/4, a larger stress in vertical 
direction. The end of the load paths in the stress space 
maps into a strain envelope response de(9) in the strain 
space. Likewise, the end of the unload paths map into 
a plastic envelope response de p {9). The yield direction 4> 
can be found from this response, as the direction in the 
stress space where the plastic response is maximal. This 
is close to 6 = 90° in this case. The flow direction ip is 
given by the direction of the maximal plastic response in 
the strain space, which is close to 70° in this example. We 
found that these directions do not agree, corresponding 
to a non- associated flow rule as it is observed in exper- 
iments on realistic soils JlTj]. If one approximates the 



plastic envelope response by its projection in the flow di- 
rection, the elastic response can be written as the simple 
form 



de p 



da) 



(14) 



x 10 4 






A 




\ \ 

- 


' / ,•••" 


V" / 


\\ y 
V' 





dp/p 



xlO 



8 








90 

S7.5 










112.5 






6 








1 / 45 


4 




135/ 






22.5 


2 


157i 










J- 

V 






\ 1 //// 




I 





180; 




/ \ 




37.5 


-2 


202.5. 






\ 315 




-4 


!2S*'' 


















292.5 








247.5 


270 






-4 -3 


-2 


V 


1 2 3 


4 

K10"' 



FIG. 2. Stress- and strain-relations resulting from the 
load-unload test, grey lines represent the paths in the stress 
and strain spaces. The dotted line gives the strain envelope 
response and the solid line is the plastic envelope response. 

The so-called hardening modulus h is the ratio between 
the modulus of the maximal plastic strain and the modu- 
lus of the incremental stress. The function (x) is defined 
as zero if x < 0; Otherwise it is valued to x. The unit 
vectors ip and 4> define the flow direction and the yield di- 
rection, respectively. The vector da defines the direction 
and magnitude of applied load. 
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FIG. 3. Elastic response de e and plastic response de p for different levels of pressure p and shear stress q. The envelope 
responses result from the application of different loading modes with \da\ — 10 _4 p. The elastic response, calculated from Eq. 

has a centered ellipse as envelope response. The plastic strain response lies almost on a straight line oriented in the flow 
direction. This direction does not correspond with the yield direction (dashed lines) . The dash-dotted line represents the failure 
surface. When a stress value above this line is applied the system fails. The solid line represents the plastic limit surface. The 
latter is obtained connecting the points where the plastic deformation diverges. (That means de p /\da\ — > oo) The dashed line 
q — 0.18p represents the limit of the zone where the isotropic elasticity assumption is valid. The plot shows the yield and the 
flow directions as a function of the deviatoric stress ratio q/p, evaluated from the average over five different samples. These 
calculations are made with fixed k n — 160 MPa. 



Both elastic and plastic envelope responses are calcu- 
lated from different stress values. The results are shown 
m Fig. | The elastic response, calculated from Eq. (I13J) , 
has a centered ellipse as envelope response for all the 
cases. Since the direction of this response does not al- 
ways correspond to the direction of the volumetric strain 
increment, the general form for the elastic stiffness must 
be written as 



de e = 



1 — v —a 
—a 1 + v 



da. 



(15) 



Here E and v are the classical parameters of the elastic- 
ity, i.e. the Young modulus and the Poisson ratio. They 
are not material parameters because they depend on the 
stress state we take. Moreover, an additional variable a 
must be included in this relation, taking into account the 
anisotropy of the elastic response. A limit of isotropy is 
found around q = 0.18p (See Fig. ||). Below of this line 



the parameters of elasticity are almost constants with 
a w 0. Above this limit the stiffness decreases as a re- 
sult of open contacts, giving rise to an anisotropic elastic 
response. 

In a previous work [ ^3j , the elasto-plastic quantities 
resulting form Eq. ( p^ ) and Eq. ( [l5| ) have been evaluated 
as a function of the stress state. Since the mechanical re- 
sponse of soils depends not only on the initial stress state 
but also on the way how this state is reached |L(|, these 
results are only valid in the case of monotonic load. In 
the classical theory of elasto-plasticity, the dependence of 
the strain response on the history of the deformation is 
described by the evolution of the so-called yield surface. 
This surface encloses a hypothetical region in the stress 
space where only elastic deformations are possible |l8|] . 

We attempt to detect the yield surface by using a stan- 
dard procedure proposed in experiments with sand |t7l . 
Fig. [| shows this procedure: initially the sample is sub- 
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ject to isotropic pressure. Then the sample is loaded in 
the axial direction until it reaches a the yield-stress state 
with pressure p and shear stress q. Since plastic deforma- 
tion is found in this stress value, the point (p, q) can be 
considered as a classical yield point. Then, the classical 
theory assumes the existence of a yield surface contain- 
ing this point. In order to explore the yield surface, the 
sample is unloaded in the axial direction until it reaches 
the stress point with pressure p — bp and shear stress 
q — 5p inside the elastic regime. Then the yield surface 
is constructed by taking different directions in the stress 
space for re-loading. In each direction, the new yield 
point must be detected by a sharp change in the slope in 
the stress-strain curve, indicating plastic deformations. 



, x 10 



yield surface 
Start of yielding 



yield point 




FIG. 4. Experimental procedure used to obtain the yield 
surface. Load-unload-reload tests are performed, and the 
points in the reload path, where the yielding begins are 
marked. The yield function is constructed by connecting these 
points. 



Fig. H shows the strain response taking different load 
directions in the same sample. If the direction of the 
reload path is the same as that one of the original load 
(45°), we observe a sharp decrease of stiffness when the 
load point reaches the initial yield point, which corre- 
sponds to the origin in Fig. ||. However, if we take a dif- 
ferent direction of re-loading, we find that the decrease 
of the stiffness with the loading becomes smooth. Since 
there is no straightforward way to identify those points 
where the yielding begins, the yield function, as it was 
introduced by Drucker & Prager [ fl8| in order to describe 
a sharp transition between the elastic and plastic regions 
is not consistent with our results. 




FIG. 5. Strain response resulting from the preparation of 
the sample according to Fig. ^. The solid lines show the 
strain response for different loading directions. The dashed 
contours connect the strain increments values with the same 
value of | Act | . The Figure shows that in any reload path dif- 
ferent to that one of the original load the yielding develops 
continuously. Thus, it is not possible to distinguish an elastic 
regime. 



IV. CONCLUDING REMARKS 

The incremental elasto-plastic response of a Voronoi 
tessellated sample of polygons has been examined in the 
framework of the classical theory. The resulting consti- 
tutive relation leads to non-linear, anisotropic elasticity, 
where the classical parameters of elasticity, the Young 
modulus and the Poisson ratio, are not material con- 
stants. The plastic response reflects the non-associated 
features of realistic soils. Here the classical analysis of 
Drucker & Prager is not applicable, because it is not 
possible to determine an elastic regime. 

Future work will be oriented to a micro-mechanically 
based description of these elasto-plastic features. Since 
the mechanical response of the granular sample is rep- 
resented as a collective response of all the contacts, it 
is expected that the macro-mechanical response can be 
completely characterized by the inclusion of some field- 
variables, which contain the information about the micro- 
structural arrangements between the grains. Some statis- 
tical variables like the fabric tensor have been included 
as internal variables pIHI]- This description, however, 
does not seem to offer a complete characterization of the 
constitutive response. More salient aspects, such as the 
yielding of the contacts, and the fluctuations of the stress 
inside the granular material, might offer a more complete 
set of internal variables for the description of the macro- 
scopic state. This is, in our point of view, an important 
challenge in the future. 
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